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Abstract 

Repeated pesticide contaminations of lentic freshwater systems located within agricultural landscapes may affect 
population evolution in non-target organisms, especially in species with a fully aquatic life cycle and low dispersal ability. 
The issue of evolutionary impact of pollutants is therefore conceptually important for ecotoxicologists. The impact of 
historical exposure to pesticides on genetic divergence was investigated in the freshwater gastropod Lymnaea stagnalis, 
using a set of 14 populations from contrasted environments in terms of pesticide and other anthropogenic pressures. The 
hypothesis of population adaptive divergence was tested on 1 1 life-history traits, using Q ST -F ST comparisons. Despite strong 
neutral differentiation (mean F ST = 0.291), five adult traits or parameters were found to be under divergent selection. 
Conversely, two early expressed traits showed a pattern consistent with uniform selection or trait canalization, and four 
adult traits appeared to evolve neutrally. Divergent selection patterns were mostly consistent with a habitat effect, 
opposing pond to ditch and channel populations. Comparatively, pesticide and other human pressures had little 
correspondence with evolutionary patterns, despite hatching rate impairment associated with global anthropogenic 
pressure. Globally, analyses revealed high genetic variation both at neutral markers and fitness-related traits in a species 
used as model in ecotoxicology, providing empirical support for the need to account for genetic and evolutionary 
components of population response in ecological risk assessment. 
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Introduction 

Understanding the causes of genetic divergence within species is 
a central goal in evolutionary biology. Population evolution is also 
a question of importance in environmental sciences and conser- 
vation biology, and, as a consequence, evolutionary principles are 
increasingly taken into account in the management of biodiversity 
and ecosystems [1-3]. Indeed, evolutionary change can be rapid, 
especially in the current environmental context. Besides evolu- 
tionary textbook cases of insecticide resistance [4] , environmental 
stress may trigger fast evolutionary responses [5], and rapid 
adaptive responses have been documented in ecological situations 
involving anthropogenic disturbance [6,7]. Furthermore, it is 
noteworthy that the rate of phenotypic change can itself increase 
in human-disturbed and toxic environments [8,9]. In spite of this, 
the potential evolutionary impact of chemicals and toxic 
substances released by human activities into the environment is 
still not considered in current procedures of ecological risk 
assessment. However, strong arguments for mainstreaming such 
impacts have been formulated for a long time [10-14]. Beyond 
mutagenic compounds affecting the germline, pollutants have 
proved to be potential sources of evolutionary impact, as 



supported empirically by an ever growing number of studies 
[15-23]. 

Agriculture is probably one of the main causes of human- 
induced microevolution by selection, as a direct consequence of 
genetic resource management, e.g., domestication, selection 
programs on productive species, or side effects of intensive 
agricultural practices [24-26]. Furthermore, in the case of 
pesticides, as substances rarely reach their target species exclu- 
sively, non-intentional evolutionary effects are to be expected. 
Arguably, a substantial loss of biodiversity has resulted from the 
intensification of arable agriculture over the last five decades [27]. 
Agriculture, combined with industrial and domestic activities, uses 
more than one-third of the Earth's accessible renewable freshwater 
(i.e., approximately 4,430 kmVyear in 2006) [28] and often leads 
to water contamination. About 140 million tons of fertilizers and 
several million tons of pesticides are applied each year [29]. 
Natural populations are thus exposed to a diversity of pressures 
that vary both spatially and temporally, and a succession of wildlife 
mortality events has marked the evolution of industry and 
agrichemical practices since the early 20 th century [30]. 

With respect to human-induced stress, some conditions are 
expected to increase the risk of genetic change in natural 
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populations. For instance, freshwater lentic habitats located within 
agricultural landscapes are repeatedly contaminated by complex 
pesticide mixtures, through various modes of transfer from the 
treated parcels (including aerial drift, run-off, and drainage; see 
[31]). Although a given pesticide or its metabolites may be only 
transiently present in water, the succession and combination of 
different molecules are expected to lead to high environmental 
variation and recurrent stress. At the population level, if 
demographic stochasticity is triggered by stressful episodes, 
random genetic drift and inbreeding will increase, and impede 
local adaptive processes [32]. Indeed, selection efficiency is 
positively correlated to effective population size [33], and 
population response to selection is weaker under high inbreeding, 
due to linkage disequilibrium and selection interference among 
loci [34-36]. Also, random drift load will accumulate locally and 
increase the risk of population extinction [33,37]. Furthermore, 
population adaptive potential is also predicted to decrease in 
stressful and changing environments, especially due to genetic 
erosion [38]. These processes will be exacerbated under limited 
gene flow, which may result from low dispersal abilities (fully 
aquatic organisms, e.g., molluscs, crustaceans) or opportunities 
(weak connectivity among occupied sites, e.g., ponds, temporary 
ditches) [39]. Also, population tolerance to environmental stress 
and change will decrease as inbreeding increases, given that 
inbreeding depression generally worsens under stress [40]. 
Interestingly, environmental heterogeneity in selection also 
increases inbreeding depression [41]. 

Alternatively, faced with stress, populations may adapt in two 
ways. First, they may adjust through phenotypic plasticity, which 
should in principle have no evolutionary consequences, although 
plasticity may itself be under selection [42,43]. Second, genetic 
adaptation may occur, provided that adapted alleles are present 
and selection intensity is strong enough to overcome stochastic 
interference [44]. 

To sum up, human-induced environmental stress and hetero- 
geneity may alter population evolutionary trajectory in varied and 
complex ways, the study of which is anything but an easy task. 
Processes described above will also have consequences on 
population genetic divergence: both random genetic drift and 
local adaptation will cause population divergence, whereas 
environmental stress and heterogeneity may either act in the 
same way (random drift or varying selection regimes) or on the 
contrary, trigger similar adaptive responses among populations 
(uniform selection acting on general stress responsive pathways or 
traits). 

One method commonly used to decipher population genetic 
divergence in terms of selective and neutral processes, is the Qst- 
F ST approach [45,46]. The general principle is to partition the 
total genetic variance at neutral markers and in quantitative traits 
into within- and between-population components, using Wright's 
F ST index, and its quantitative analogue Qk,t> respectively. Under 
a purely additive determinism, if quantitative traits evolve 
neutrally, they should lead to similar levels of differentiation as 
those inferred from neutral markers (Qst = F st). Basically, (5jj T 
and F S t are compared to test the null hypothesis of neutrality. 
Differences between the two indices indicate selection on the traits, 
with two possibilities: homogenizing (uniform stabilizing) selection 
towards a unique fitness optimum is revealed by ^st- > !2st> 
whereas divergent selection (occurrence of different fitness optima 
among populations) is reflected by the reverse relationship. This 
principle has been widely used in the last decade to address 
population adaptive divergence (see [46] and references therein). It 
has notably revealed population adaptive divergence as a function 



of habitat [47], environmental acidification [48], or soil contam- 
ination by zinc [49]. 

In the present study, we applied the Qst-^st method to test the 
hypothesis that anthropogenic pressure may act as a selective force 
and lead to population adaptive divergence, using a set of 14 
natural populations of the freshwater snail Lymnaea stagnalis 
(Mollusca, Gastropoda, Panpulmonata, Hygrophila; previously 
classified in the Sub-Order Basommatophora). Due to its habitat 
characteristics (lentic systems, often close to agricultural zones) and 
its limited dispersal ability, this species is a good model to test non- 
intentional evolutionary effects of pesticides. Moreover, L. 
stagnalis belongs to an ecologically important group, which plays 
a major role in the transfer of energy across food webs, and can 
represent up to 20-60% of the total biomass of macro- 
invertebrates in some freshwater ecosystems [50]. Pesticides and 
other identified anthropogenic pressures were described qualita- 
tively and considered as contributing to a global and potentially 
toxic pressure on non-target organisms. In this context, directional 
selection associated with a particular pesticide was considered 
irrelevant, and thus not specifically focused. Conversely, we 
hypothesized that global chronic exposure to various cocktails of 
pesticides may induce selective response to temporary but 
recurrent stress. The study aimed at addressing the following 
questions: (i) is agricultural (pesticide) pressure able to trigger 
adaptive processes in non-target organisms? (ii) what is the relative 
strength of such a pressure, compared to other environmental 
components, including natural and other anthropogenic factors? 
Indeed, as natural populations in a given level of pesticide pressure 
are not eco-evolutionary replicates of each other, any other 
relevant environmental factor should be also considered. Under 
this rationale, we characterized environments in terms of 
"pesticide pressure", "other anthropogenic pressures" (roads, 
urbanized zones), "global pressure" (combination of agricultural 
and non agricultural anthropogenic pressures), and "habitat" 
(pond, ditch, channel, as physical features expected to affect 
population size and isolation). Population divergence patterns 
were investigated for a set of life-history traits and compared to 
neutral genetic differentiation. 

Materials and Methods 

Population sampling and characterization 

Adult individuals were collected from 14 natural populations 
during reproduction, in the course of a single campaign (5-12 July 
201 1) (Fig. 1). Sampling was conducted in public locations that did 
not require specific authorization, and did not involve endangered 
or protected species. Locality information and sampling charac- 
teristics are summarized in Table 1. Pesticide and other anthro- 
pogenic pressures were estimated from land-use patterns observed 
in the immediate surroundings of sampled sites, and using Google 
Earth. In the present study, we were not interested in an instant 
environmental concentration of pesticides at the time of sampling. 
On the contrary, under the tested hypothesis, we assumed 
historical contamination (repeated or chronic) as hypothetical 
selective force driving phenotypic evolution. Our strategy built on 
the use of landscape geographic information and agricultural land- 
use data. These have proved to be relevant tools to predict water 
body contamination by pesticides [51] and to assess pesticide 
exposure in the field (ditches [52]; ponds [18,53,54]). Moreover, in 
such aquatic systems, most variation may be explained by close 
land-use patterns (within a 100 m radius area), as demonstrated 
for water quality and vegetation complexity [55] . Consistently, we 
characterized each population's environment in terms of percent 
coverage by forest or moorland, pasture, crop (including potatoes, 
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corn, orchards, bulb plants), and urban zone (Table 1). By crossing- 
field observations and Google Earth satellite views with informa- 
tions obtained from farmers on their practices, we estimated 
coverage proportions using the software Image] (http://rsbweb. 
nih.gov/ij/). When available, pesticide concentration data in 
surface water near the sampled sites were found to be globally 
consistent with our classification (http://81.93.58.66/ 
bma_nieuw/begin.html; http:/ /www.milieurapport.be/). Finally, 
populations were also classified according to four criteria: habitat 
(H: pond, channel, ditch), pesticide pressure (PP: two levels, low vs 
high), other anthropogenic pressure (OAP: two levels, low vs high), 
and global environmental pressure, as a combination of pesticide 
and other anthropogenic pressures (GEP: three levels). 

Molecular analyses 

DNA was chelex extracted from haemolymph or foot tissue 
from 399 wild-caught adults (14 to 33 snails per population). 
Neutral genetic variation was assessed at 12 microsatellite loci, i.e., 
A2, A112, Bl 17 [56], 2kll and 2k27 (Genbank accession: 
EF208747-EF208748 [57]), and EMLS04, EMLS13, EMLS21, 
EMLS26, EMLS29, EMLS41, EMLS45 [58], following the 
protocol described in Besnard et al. [58]. Only individuals with 
less than three missing genotypes were retained for population 
genetics analysis. 

Population neutral genetic structure 

Mean allele number (N), allelic richness (A R ), expected 
heterozygosity H E , and observed heterozygosity H 0 , were 
calculated with Genetix 4.05.2 [59]. The distribution of neutral 
genetic diversity within and among populations was estimated 
from Weir and Cockerham's estimators of Wright's F indices [60] 
using Fstat 2.9.3.2 [61]. Departures from HWE (heterozygote 
excess or deficiency) and linkage disequilibria were tested using 
Genepop 4.0.10 [62]. Population differentiation was tested with a 
permutation test, in which genotypes were permuted among 
samples (not assuming HWE within samples; see [61]). As L. 
stagnalis is a self-fertile hermaphroditic organism, the selling rate 
was estimated per population and statistically compared to zero 




Figure 1. Location of 14 L stagnalis populations involved in <? ST - 
F ST comparisons. Populations are coded from 1 to 14. Colours 
indicate environmental categories (white, grey, and black, for GEP, 
GEP1, and GEP2, respectively). See Table 1 for details. 
doi:10.1371/journal.pone.0106670.g001 



using Rmes [63]. Effective population size was estimated using the 
sibship assignment method, as implemented in the software 
Colony 2.0.3.0 [64] and assuming inbreeding, male and female 
polygamic mating systems, and monoecy. 

To estimate the number of genetic clusters in our dataset 
without taking into account any predefined population, we used 
Structure 2.2 [65]. Analyses were performed assuming an 
admixture model and a number of genetic clusters (k) from 1 to 15. 
Each run started with a burn-in period of 50 000 steps followed by 
300 000 Markov Chain Monte Carlo (MCMC) replicates. The 
most likely number of clusters was determined using the Ak 
statistic [66] using Structure Harvester [67]. We used 
Distruct to plot Structure output data [68]. 

The effect of environmental factors PP, OAP, GEP, H, and 
genetic cluster was assessed on population genetics parameters 
(Ap L , He, Fis, Fst) with a permutation test using Fstat. 

Common garden experiment 

Wild-caught adults were brought to the laboratory and reared 
under standard conditions at the INRA Experimental Unit U3E 
(Rennes, France), as previously described [69] . Snails were isolated 
in plastic vessels filled with 1 L de-chlorinated and charcoal filtered 
tapwater. They were fed weekly with 1 .5 g of organic salad, at 
each water renewal. Room temperature was maintained at 
20± 1°C and the photoperiod was 16L/8D. Like other basomma- 
tophorans, L. stagnalis lays eggs embedded in mucous enveloppes 
which are deposited and fixed on available substrates (sediment, 
vegetation, tank walls, etc.). For a given snail, reproduction was 
followed during 1 4 days after the first clutch laid in the laboratory. 
A total of 228 snails over 399 reproduced after 3 months (57%). 
Families were characterized at various life-history traits (individual 
growth, female reproduction, hatching success), measured on the 
laboratory-born progeny (Gi) as illustrated on Fig. 2. G t snails 
were reared at 20±1°C, under a 14L/10D photoperiod. 

Gf rearing conditions 

From hatching to the age of 1 1 9 days (roughly corresponding to 
the end of juvenile growth), Gi snails were reared as groups 
reflecting clutch origin (n = 427; 1 to 2 clutches per maternal 
parent). As individuals grew, rearing conditions were adjusted in 
terms of water volume, individual density, and food supply (see 
Table SI). Vessels were regularly moved in a randomized way. 
From the age of 1 1 9 days, 7 individuals per clutch were randomly 
chosen and marked with a honey bee mark (model FC075; 
diameter: 2.3 mm, weight: 1.8 mg; Ickowicz Apiculture, Bollene, 
France), and further reared in 30 L tanks (80 individuals per tank), 
in which they were fed twice a week with organic salad (0.5 g per 
snail), and from which 2/3 of water was renewed weekly. 

Life history traits: individual growth 

From hatching to the age of 56 days, individual size (as inferred 
from shell height) was measured every two weeks on four 
randomly chosen juveniles per group. From the age of 63 days, 
all reared snails were individually measured. Measurements were 
performed with a stereomicroscope fitted with an ocular micro- 
meter until day 63 , and a digital calliper afterwards. From the age 
of 1 19 days, three G! snails per family (n = 492) were individually 
followed for growth and reproduction. For these 492 snails, size 
was measured at the age of 1 77.7 1 ±5.41 days, 191.71 ±5.41 days, 
205.7±5.40 days, and 226.5 1 ±5.56 days (owing to variation in 
clutch age, and because measurements were performed at fixed 
dates after the age of 1 19 days), leading to a total of 12 values per 
individual. Growth was modelled using the Gompertz's model: 
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where L t is the shell height at time t, A the asymptotic shell length, 
b a scaling factor related to shell height at t = 0, and k reflects the 
growth rate. Using a preliminary subset of 20 individuals, this 
model performed better than the von Bertalanffy model (lower 
AIC), and was thus preferred for further growth analysis. Growth 
parameters were compared among families and populations, as 
well as size at hatching and size at 1 19 days. 

Life history traits: reproduction 

Once clutches started to be observed in all aquaria, the 492 
snails followed for growth were isolated in 200 mL plastic vessels 
and fed weekly with 1 g of organic salad immediately after water 
renewal. Several female reproductive traits were measured: time to 
first oviposition under isolation (with a censoring limit of 30 days), 
ability to lay eggs (as the proportion of reproductive snails per 
family), number of clutches and eggs laid during two weeks, clutch 
size (number of eggs per clutch), and clutch hatching rate. 

Statistical analyses 

All traits were analyzed with generalized linear mixed effect 
models (GLMM, R-package lme4) [70], with appropriate error 
distribution (Gaussian for normal data, Poisson for count data, and 
binomial for proportions). When necessary, data were log- or 
BoxCox-transformed, and covariates were included in the model 
(see Table S2). The model structure was: Y ~ factor 1 + factor! + 
factor3 + covariatel + covariate2 +...+ (1 {population/family) . 

Fixed effects were tested by model comparison using a log- 
likelihood ratio test. Fixed factors were: genetic cluster (see 
Structure analysis), habitat, and environmental pressure (either 
pesticide, other anthropogenic, or global pressure). Family was 
nested within population, and both were treated as random 
factors. Global tests were followed by post-hoc pairwise compar- 
ison tests (Tukey and non parametric equivalent). All statistical 
analyses were performed using R 2.14.0 (R Core Team 2012). 

Qst -F ST analysis 

Qst was computed for each trait using the equation [71;72]: 



(l+f) V b + 2 V w 



Where / is the inbreeding coefficient, Vj, the between population 
genetic variance and V w the within population genetic variance. V\, 
was obtained directly from the population variance component, 
while for V w , as only the maternal origin could be assessed, family 
variance was used under the full-sib design (2VJ) [73]. Within- and 
between-population variance components were estimated for each 
of the 1 1 studied traits using restricted maximum likelihood 
(REML) under a simplified mixed model (Y~l+ (1 \ pop/fam)). 

Population quantitative divergence was tested using the method 
of Whitlock and Guillaume [74], which accounts for variation 
among neutral loci used to estimate Fst, and for sampling error 
and variation in evolutionary history in Qst estimation. The 
method has lower type 1 error rate and better statistical power 
than previous tests, and seems particularly suited when neutral 
differentiation is high [46], as in L. stagnalis [58,75]. Briefly, a 
bootstrapping procedure is applied to compare observed values of 
the difference Qjt - ^st to the (5st - ^st distribution expected 
under the neutral hypothesis, as derived from F ST distribution and 
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Figure 2. Schematic overview of the common garden experiment used to investigate population divergence in L. stagnalis. 

doi:1 0.1 371 /journal.pone.01 06670.g002 



variance components [60] and Qst neutral distribution predicted 
from the Lewontin-Krakauer distribution. Observed values were 
considered significant when they fell outside the 95% confidence 
interval of the neutral distribution [74]. 

Results 

Population genetics 

On 399 individual genotypes, 44 were discarded because of 
missing or unreadable data at more than three loci. Overall, the 
observed number of alleles per locus varied from three to 29, mean 
allelic richness ranged from 1.1 to 9.6 (Table S3), and genetic 
diversity per locus and sample varied from 0 to 0.86. Significant 
genotypic linkage disequilibrium was observed in 40 out of 656 
comparisons. However, after Bonferroni correction, none of these 
remained significant. Fis-values varied widely across loci and 
samples, ranging from -0.238 (EMST41 in BUX) to 1.000 (2k27 
and Bl 17 in DET; EMLS13 in HED; EMLS21 in BAA) (Table 
S3). Genetic parameters estimated per sample over loci reflected 
discrepancies between populations (Table 1). All populations 
except BUX and OUD were found significantly inbred (mean 
F IS -value =0.219, 95% CI [0.147; 0.287]). The selfing rate was 
significandy different from zero in three populations (OUD, DET, 
and AGA). Effective population size Ne was significantly larger in 
sites exposed to high pesticide pressure (PP 1 > PP0; Kruskal- W allis 
test, P = 0.007). 

Population differentiation was generally high (mean F^j-value = 
0.291, 95% CI [0.248; 0.334]), as also reflected by pairwise 
estimates, which were all significant (Table S4), except between the 
two geographically close populations KUI and EMM. Interestingly, 
the four pond populations had the longest branches on the NJ-tree, 



indicating stronger differentiation associated with this habitat 
(Fig. 3A). From Structure analysis, individual genotypes clustered 
into two groups with the highest likelihood, and this was confirmed 
by the Ak statistic (Fig. 3B). Clusters corresponded to two 
geographic regions (10 western vs. four eastern populations). 
Therefore, the effect of "genetic cluster" was tested on life-history 
traits. No difference in level of genetic diversity, population 
inbreeding or population differentiation was found between the 
two clusters (Table 2). 

Population genetic parameters were affected by most environ- 
mental factors (Table 2). Genetic diversity tended to be greater 
under higher environmental pressure (GEP2, PP1, OAP1), 
although the difference was significant under pesticide pressure 
only {A K , P = 0.007; H E , P = 0.008). Population differentiation was 
significandy affected by global environmental pressure (GEP, 
P = 0.016) with a lower value among populations exposed to the 
highest level of pressure. With regard to specific pressures, lower 
differentiation was associated with pesticide as well as with other 
anthropogenic pressures (PP, P = 0.001; OAP, P = 0.031). Habitat 
affected all tested parameters, when ponds were compared to ditch 
and channel populations: pond populations were significantly less 
variable {A R , P = 0.025; H E , P= 0.038), more inbred (F IS , 
P = 0.034) and much more differentiated than were other 
populations (F ST(l , ond) = 0.502 vs. F ST (ditch + channel) = 0.227, 
P = 0.002). 

Life history variation 

Most G[ traits showed significant heterogeneity at both 
population and family levels (Tables S2 and S5). No effect of 
genetic cluster was detected, excepted on early size, which was 
larger in eastern populations (size at hatching, P = 0.019; growth 
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parameter b, P = 0.005; Table 3). Among environmental factors, 
habitat was the most effective, with six traits or parameters 
significantly different between systems. Growth parameters 
indicated that body size at t 0 increased significantly from ponds 
to channels and from channels to ditches, whereas asymptotic size 
was larger in ponds than in channels and ditches (Fig. SI). With 
respect to reproduction, snails from ponds had a significantly lower 
reproductive activity than those from ditches and channels: lower 
ability to lay eggs, and in reproductive snails, longer time to 
oviposition. Clutch size decreased significandy from ponds to 
ditches and from ditches to channels, and clutch hatching rate was 
lower in pond and ditch populations than in channel ones. 

Compared to aquatic system, human pressure appeared much 
less effective. First, pesticide pressure had no effect, except a 
marginal one on ability to lay eggs, found greater in exposed 
populations (P = 0.069). Second, other anthropogenic pressures 
affected growth parameter b (significandy larger size at t 0 under 
OAP1) and hatching rate (significantly lower under OAP1) (Fig. 
S2). Interestingly, lower hatching rates were already observed on 
the clutches laid by G n snails inhabiting OAP1 sites (720 clutches, 
GLMM test: x 2 = 1 1-861, P<0.001). Third, global pressure 
affected b (significandy larger snails at to under GEP2), and 
hatching rate, which was significantly impaired in exposed 
populations (GEP1 and GEP2) relative to reference ones (GEP0). 

Several life history traits correlated with population genetic 
characteristics. First, fecundity was positively correlated with H E 
(r = 0.524, P = 0.03 1 ; Fig. 4A), and with A R (r = 0. 489, P = 0.044). 
Second, a negative correlation was found between population 
inbreeding and G, fecundity (r= -0.845, P<0.001; Fig. 4B), 
ability to reproduce (r= —0.664, P = 0.006), and clutch hatching 
rate (r= -0.634, P = 0.009) (see Table 1 and Table S5). 

Qst-Fst analysis 

Results are summarized on Table 4. Over the set of 1 1 studied 
traits or parameters, four traits were found to evolve consistently 
with neutral expectations (hatching size, growth parameter k, 
number of eggs, clutch size), whereas homogenizing selection was 
indicated for two traits (growth parameter b, clutch hatching rate), 
and divergent selection for five traits (size at 1 19 days, asymptotic 
size A, ability to lay eggs, time to oviposition, number of clutches 
laid). 

Discussion 

The hypothesis that human activities may induce population 
adaptive divergence in L. stagnalis was addressed under a 
common garden experiment involving life history traits and a set 
of natural populations from contrasted environments. Five traits or 
parameters over 1 1 were found under divergent selection, despite 
a strong neutral genetic structure. The occurrence of two genetic 
clusters (Eastern and Western populations) had only a weak 
influence on quantitative genetic divergence, suggesting that 
population level was the adequate scale to draw selection 
inferences [76]. Only size at hatching and its corresponding 
growth parameter (b) were affected by this factor (larger hatchlings 
to the East). Given that early size was found under uniform 
selection, it was assumed that the genetic divergence associated to 
geography did not interfere with selection putatively associated 
with local environmental pressures. 

Selection versus random genetic drift 

All traits found to diverge adaptively were adult traits, related to 
late expressed growth (sub-adult size, asymptotic size A) and to 
reproduction (ability to reproduce, time to oviposition after 




Figure 3. Population genetic differentiation in L. stagnalis. (A) 

Bayesian assignment probabilities in Structure analysis, for k = 2 clusters. 
Each bar represents an individual. Bar colour indicates the posterior 
probability that the individual belongs to the cluster of that color. (B) 
Unrooted Neighbour-Joining tree based on population pairwise F ST 
values (see Table S4). 
doi:1 0.1 371/journal.pone.01 06670.g003 

isolation, number of clutches). Conversely, uniform selection was 
the most likely evolutionary hypothesis for early survival and size 
(hatching success and growth parameter b, which is related to size 
at t 0 ). Globally, inferred patterns, i.e., early (late) traits under 
homogenizing (divergent) selection, are in lines with those found in 
another freshwater snail [47,77]. 

Traditional statistical analyses performed on the Gi revealed 
that, to the exception of a marginal effect of pesticide pressure on 
ability to lay eggs (P = 0.069), habitat was the only factor affecting 
most traits found under divergent selection. Indeed, pond 
populations were characterized by larger ultimate size (growth 
parameter A), lower ability to lay eggs, and longer time to 
oviposition after isolation, compared to ditch and channel ones. 

Congruent with L. stagnalis low dispersion and weak habitat 
connectivity, pond populations were particularly differentiated 
(especially CAS and DET, pairwise F^t) and presented generally 
lower genetic diversity and higher inbreeding than channel and 
ditch populations. Moreover, compared to channel populations, 
pond and ditch populations were also more inbred (P = 0.009). 
These aquatic systems are likely to imply stronger population 
isolation, even among ditches (due to frequent drought). A similar 
influence of habitat was emphasized in another freshwater snail, 
Physa acuta [78]. Thus, in L. stagnalis, although divergent 
selection was detected for traits also affected by habitat type, 
genetic drift is likely to highly contribute to divergence among 
habitats. This result is in line with the detection of significant drift 
load in experimental populations occupying outdoor close 
mesocosms [69]. 

In the common frog Rana temporaria, habitat fragmentation 
correlates with low genetic diversity and high differentiation, and 
negatively impacts tadpole body size [79], which is a critical trait 
determining individual fitness [80]. In L. stagnalis, although 
ultimate body size was greater in ponds, the relation between size 
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Table 2. Comparison of L stagnalis population genetic parameters as a function of different grouping factors. 



Population grouping 








Fst 


GEPO 


2.898 


0.367 


0.249 


0.431 


GEP1 


2.725 


0.353 


0.190 


0.450 


GEP2 


3.599 


0.468 


0.214 


0.197 


P (GEP) 


0.073 


0.094 


0.487 


0.076 


PPO 


2.824 


0.362 


0.226 


0.413 


PP1 


3.599 


0.468 


0.214 


0.197 


P(PPi 


0.007 


0.008 


0.773 


0.001 


OAPO 


2.898 


0.367 


0.249 


0.431 


OAP1 


3.337 


0.438 


0.210 


0.248 


P(OAP) 


0.219 


0.178 


0.329 


0.031 


H Pond 


2.678 


0.343 


0.282 


0.502 


H Ditch 


3.300 


0.436 


0.246 


0.221 


H Channel 


3.550 


0.449 


0.166 


0.246 


P(H) 


0.062 


0.095 


0.076 


0.001 


H Pond 


2.678 


0.343 


0.282 


0.502 


H Channel+Ditch 


3.425 


0.443 


0.203 


0.227 


P (P vs D+C) 


0.025 


0.038 


0.034 


0.002 


H Pond+Ditch 


3.024 


0.398 


0.258 


0.336 


H Channel 


3.550 


0.449 


0.166 


0.246 


P (P+D vs Q 


0.116 


0.306 


0.009 


0.271 


H Pond+Channel 


3.162 


0.409 


0.203 


0.330 


H Ditch 


3.300 


0.436 


0.246 


0.221 


P (P+C vs D) 


0.686 


0.641 


0.285 


0./78 


GC West 


3.187 


0.422 


0.229 


0.244 


GC East 


3.272 


0.410 


0.196 


0.162 


P(GC) 


0.827 


0.837 


0.426 


0.393 


A R is allelic richness, H E expected heterozygosity, F !s 


inbreeding coefficient, and F ST 


differentiation index. 







GEP holds for Global environmental pressure, PP for Pesticide pressure, OAP for Other anthropogenic pressure, H for Habitat (P= Pond, D= Ditch, C= Channel}, GC for 
Genetic cluster. Permutation-based statistical test (P value, 1000 permutations) 
doi:1 0.1 371 /journal.pone.01 06670.t002 

and fitness is not necessarily positive, because of possible energetic 
trade-offs between growth and reproduction [81]. Indeed, two 
pond populations (CAS and DET) exhibited both highest mean 
asymptotic size and lowest fecundity (see Table S5). 

The Q$t~Fst based test concluded to neutral divergence for size 
at hatching, growth rate (parameter k), and two reproductive traits, 
clutch size and early fecundity (number of eggs laid during 1 4 days 
following first clutch since isolation). On the latter, it seems 
surprising to find no hint of selection, either uniform or divergent, 
although this may simply reflect low heritability of fitness traits 
[82]. In support of this, it is worth noting that fecundity is 
extremely high in L. stagnalis (from adulthood to death, i.e., about 
one year under standard conditions, a snail produces more than 1 0 
progeny per day), suggesting the possibility of high variation 
without drastic fitness consequences. Consistently, population 
growth rate was found to be insensitive to fecundity under 
experimental conditions (Leslie matrix modelling [83]). Alterna- 
tively, apparent neutral evolution might result from a bias in Q,st 
estimation. Under genetic drift, dominance and epistasis bias Qst 
estimates downward [84—86]. As Qst is traditionally estimated 
without accounting for these non-additive sources of variation, and 
given the strong observed global F ST value, divergent selection 
might have been underestimated for fecundity. On the other side, 



a hypothetical upward bias in £)st estimation might result from 
too high a mutation rate at the markers used to infer Fst> as 
compared to the mutation rate of loci encoding the studied 
quantitative trait [87,88]. This hypothesis is invalidated by 
observed Fst values, which were rather high and significant, 
although they may also suffer from a downward bias. Finally, gene 
flow may well lead to an apparent pattern of adaptive divergence 
at a trait under uniform selection, if this trait is correlated to a trait 
under divergent selection [89]. Again, high population neutral 
differentiation makes this hypothesis unlikely. Besides methodo- 
logical causes of bias, it might be also mentioned that (1) true 
lifespan fecundity was not measured, and (2) as a source of 
uncontrolled variation, no distinction could be made between 
outcrossing and potential selfing in the Gi. 

As reported above, homogenizing selection was indicated for 
early traits, which is expected to reflect the occurrence of a 
uniform fitness optimum between the studied populations. <2st 
estimation was based on the Gi of wild-caught adult snails, with 
the exception of survival at hatching (measured on G2). Gi families 
correspond hence to the hermaphroditic equivalent of isofemale 
lines, as used for broad-sense heritability estimation [73]. 
Therefore, these "broad-sense" <2st estimates are potentially 
biased by genetic non additive sources of variation (dominance, 
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Figure 4. Correlation plot between 14-days fecundity and (A) population expected heterozygosity, (B) population inbreeding, as 
based on 14 L stagnalis populations. 

doi:10.1371/journal.pone.0106670.g004 



epistasis), and by maternal effects. The latter are likely to operate 
early in life and tend to recede with age [90] . Hence, in the present 
study, traits found under uniform selection (early expressed traits) 
might be more prone to such a bias than traits found under 
divergent selection (late expressed traits). However, as complex 
traits seem to vary mostly additively [91,92], we hypothesized, as 
done in most Qst-^st comparisons [46], that neglecting these 
sources of variation would have little consequences. 

Alternatively to the uniform selection hypothesis, the pattern 
exhibited by early traits would result from trait canalization, as a 
phylogenetically inherited characteristics [86]. Since the pattern 
was similar to that observed in Galba truncatula [77], the criterion 
of cross-species consistency (supporting the canalization hypoth- 
esis) may be met in basommatophorans, but this would need to be 
confirmed in other species. Further investigation is clearly needed 
to disentangle both causes of L. stagnalis population phenotypic 
convergence. 

Genetic diversity and fitness 

The positive correlation observed between population genetic 
variability and G! fecundity, is consistent with previous findings in 
L. stagnalis [93]. Under the assumption that Gl performances 
reflect the fitness of natural populations, the present result also 
supports the hypothesis that local drift load can be strong in this 
species, as previously suggested [69]. Therefore, under high 
random genetic drift, as indicated here, it might also be asked to 
what extent Qgx is modified by the expected stochastic evolution 
of slightly deleterious alleles as compared to pure selection-based 
predictions (random drift load [37]). 

Furthermore, the negative correlation observed between pop- 
ulation inbreeding and several traits related to fecundity (ability to 
reproduce, time to oviposition, early fecundity, and hatching rate) 
is in line with the occurrence of reduced fitness in inbred 
populations, and thus of inbreeding depression in these popula- 
tions. However, as inbreeding was found significant in most 
populations and because effective population size was usually 
small, it is suggested that random genetic drift is actually the main 
cause of these apparent correlations, through reduced genetic 
diversity and drift load accumulation caused by population 
isolation and small size [32]. This hypothesis is still strengthened 
by the fact that population inbreeding correlated negatively with 
G 2 survival at hatching (i.e., irrespective of individual inbreeding 
level). Inbreeding depression estimated under selfing relative to 



random outcrossing (ID) is particularly low in L. stagnalis [69,93- 
95], despite a clear preference for outcrossing [96,97]. Thus, 
although ID was not estimated, it seems unlikely that selfing could 
be responsible for the observed pattern. In support of this, no 
correlation was observed between population selfing rate and G! 
fecundity. 

Evolutionary impact of anthropogenic pressures 

Neutral genetic variation was inflated in populations exposed to 
anthropogenic pressures, including pesticides. This was reflected in 
terms of genetic diversity, allelic richness, and effective population 
size, and through lower genetic differentiation. These results are 
opposed to those found in experimental populations of L. stagnalis 
exposed to cocktails of pesticides [98] . More generally, they are 
not consistent with the hypothesis of increased local stochasticity 
due to anthropogenic pressure (see introduction). High genetic 
diversity may result either from the maintenance of large 
population size despite stressful conditions, or from significant 
gene flow (immigration). Alternatively, undetected population 
subdivision may also be responsible for the observed patterns, as 
subdivision is expected to maintain diversity at the global scale. 
Furthermore, at the same scale, heterozygote deficiencies are also 
to be expected (Wahlund effect). Therefore, population subdivision 
into local demes might be responsible for part of the within 
population fixation observed in the dataset. As a possible impact of 
anthropogenic pressure, this potential effect should deserve specific 
attention, e.g., through the scoring of individual location at the 
time of sampling, and using assignment tests. As already 
mentioned, the effect of anthropogenic pressures on life history 
traits was much lower than that of habitat. This was particularly 
marked for pesticide pressure, which was the primary focus of the 
study. Three possible causes for this may be discussed: 1) true lack 
of pesticide effect, 2) effects underestimated due to confounding 
environmental factors, and 3) experimental conditions too benign 
to detect adaptive divergence. 

Besides the hypothesis of pesticide innocuousness to gastropod 
populations, one explanation may be that sampling did not reflect 
the effective magnitude of pesticide pressure. However, as the 
design included areas of intensive agriculture, and due to the close 
proximity of water bodies and treated agricultural parcels, this 
seems unrealistic. Alternatively, evolutionary patterns associated 
with pesticide exposure may have been masked by other 
environmental characteristics, an idea supported by the significant 
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habitat effect observed on most studied traits. Thus, although the 
number of populations and families was close to the recommended 
values [72], and despite the ability of mixed effect models to treat 
unbalanced designs [70], our experiment may lack sufficient 
statistical power, probably due to the lack of ponds within 
agricultural zones and of channels in reference sites. Meanwhile, 
the present study emphasizes the need to take habitat into account 
in Qst-Fst comparisons [79]. 

Third, the apparent lack of effect of pesticide pressure on the 
studied traits might also result from experimental conditions, 
which were too benign. For example, the trend for enhanced 
ability to reproduce observed in the G[ from polluted areas, might 
reflect a specific response of exposed populations to overcome 
early mortality induced by chemicals (hatching rate: OAP0> 
OAP1 in Gi and G 2 ), which was unexpected to express under 
common garden conditions. In the context of evolution under 
stressful environments, trait expression would be best studied 
under a gradient of stress conditions, as also recommended for 
molecular traits and pathways [15]. Stressful conditions may 
indeed facilitate the detection of divergent selection patterns with 
respect to pesticide historical exposure, since Qst nas been shown 
to depend on the experimental environment [48,99]. 

Finally, the present study also emphasized the occurrence of 
high genetic variation both at neutral markers and fitness-related 
traits in a species used as model in ecotoxicology. This finding 
provides empirical support for the need to account for genetic 
variation in toxicity testing and, as a perspective, in future 
procedures of ecological risk assessment [12]. 

Supporting Information 

Figure SI L. stagnalis individual growth (Gj snails) 
according to population habitat. Gompertz growth curves 
based on 12 dates of measurement, performed on 492 individuals. 
(TIFF) 

Figure S2 L. stagnalis hatching rate as a function of 
anthropogenic pressure (OAP, n = 1537 clutches). 

(TIFF) 
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